A model of dispersive transport across sharp interfaces between porous materials 
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Recent laboratory experiments on solute migration in composite porous columns (Berkowitz et 
al. [H) have shown an asymmetry in the solute arrival time upon reversal of the flow direction, 
which is not explained by current paradigms of transport. In this work, we propose a definition for 
the solute flux across sharp interfaces and explore the underlying microscopic particle dynamics by 
applying Monte Carlo simulation. Our results are consistent with the findings reported in [lj and 
explain the observed transport asymmetry. An interpretation of the proposed physical mechanism 
in terms of a flux rectification is also provided. The approach is quite general and can be extended 
to other situations involving transport across sharp interfaces. 
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Berkowitz et al. [l[ have recently reported experimen- 
tal results for solute transport across the sharp interface 
between two porous materials, which show that, contrary 
to current model predictions, solute arrival times depend 
on the direction of the flow crossing the sharp interface. 
Understanding the precise physics of particle transport 
across sharp interfaces is of critical importance also for 
such diverse applications as heat flow in composite me- 
dia or the pattern evolution of morphogens govern- 
ing living tissue development The aim of this work 
is to provide a consistent physical model for the flux of 
a solute across two porous materials with different dis- 
persive properties that explains the experimental results 
presented in 

Transport processes can be conveniently described in 
terms of particles stochastically traveling within a porous 
medium. Depending on the intensity of the underlying 
velocity field, transport in porous materials may assume 
a dispersive or a diffusion-dominated character, in both 
cases leading to significant spread of an initially close 
ensemble of particles. The time evolution of this spread- 
ing is governed by a simple mass-balance equation that 
relates the time derivative of the particles concentra- 
tion c(x,t) to the divergence of the particles flux j(x,t), 
dtc(x, t) = — V- j(x, t). The interaction of the tracer par- 
ticles with the fluid flow manifests itself at the macro- 
scopic scale as a complex combination of diffusion, dis- 
persion, and time-memory effects, which are ultimately 
determined by the details of the pore geometry Q . 

Transport across an interface is a deceptively simple 
problem that has long been studied 0, (| 0, H|, and 
for which several alternative models have been proposed. 
Some of these models assume equality of fluxes and con- 
centrations at the interface 0, [l(| • Other studies relax 
the assumption on equality of concentrations, but either 
neglect advection [i[ , or consider only concentration pro- 
files [III E3 • None of these models has been so far cor- 
roborated by experimental evidence of a concentration 
jump at the interface. 

To fix the ideas, let a composite porous column of 
length L = 1 consist of two half-columns of diameter 
D, filled with 'fine' (F) and 'coarse' (C) random arrange- 



ments of spheres of diameter dp and dc, respectively, 
with dp < dc- A fluid flux, Q, can be injected from 
cither end, fine or coarse, of the column. We assume 
that the aspect ratio of the column, L/D, is sufficiently 
large to justify a onc-dimcnsional treatment. NMR ex- 
periments show that, for random packings of spheres, 
porosity n, electrical tortuosity, and diffusivity do not 
depend on the spheres' diameter, d 13(. From the fluid 
incomprcssibility condition, the mean transport velocity 
v = Q I '(nirD 2 1 '4) also does not depend on d. Moreover, 
for the Pcclet numbers considered in 0, Pe= [0.3-10], 
diffusion is negligible with respect to dispersion. The 
relevant material properties that depend on the sphere 
diameter d are the permeability k ~ gP, which charac- 
terizes the fluid flow in the pore geometry, and the dis- 
persivity length a, the ratio of the effective dispersion 
coefficient D in the Taylor- Aris regime to the pore veloc- 
ity v, which characterizes the macroscopic spreading of 
the tracer plume through the microscopic pore geometry, 
a= D/v~d dEH. 

While each layer can separately be regarded as be- 
ing homogeneous, we have now a sharp interface (i.e., 
a macroscopic heterogeneity) at x = 1/2, the junc- 
tion between the layers. With reference to the fine- 
to-coarse (F— >C) direction, we assume that a(x) can 
be represented by a Heaviside step function, i.e., a(x) 
takes the values Of, ac, and a = (ap + etc)/ '2, for 
x < 1/2, x > 1/2, and x = 1/2, respectively. An analo- 
gous variation will be assumed for the the permeability 
kp < k(x) < he- Obvious modifications hold for the 
coarse-to-fine (C— >F) flow direction. 

Concerning the particles flux, it is customary to adopt 
a Fickian advection-dispersion (AD) constitutive rela- 
tionship j AD (x, t) = v [c(x, t) - a{x)d x c(x, t)} [9]. The 
AD model predicts no difference in the arrival times at 
the column outlet, as readily shown by numerical integra- 
tion. This is however inconsistent with the experimental 
evidence [l|, where the solute transported in the F^C 
direction is observed to arrive faster with respect to the 
one transported in the C— >F direction. Similarly, the 
adoption of a Fokkcr-Planck (FP) constitutive relation- 
ship j FP (x, t) = v [(1 — d x a(x))c(x, t) — a(x)d x c(x, t)] as 
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FIG. 1: Deviations from the average Stokes velocity, 
u = u— < u >, at the interface between two-dimensional 
disks (arbitrary units). The histogram of u for the F^C 
flow direction exhibits a distinct positive skewness. Negative 
skewness is observed for the C— >F direction (not shown). 
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FIG. 2: Skew-normal jump-length pdfs p(s\x) for the case of 
F-^C flow direction. The transition region is delimited by 
two vertical dashed lines. 



in [14] predicts faster solute arrival times are for the C— >F 
flow direction and thus does not explain the results in [l[ . 

Bearing in mind these results, in the following we pro- 
pose a model for the microscopic dynamics of a pas- 
sive tracer across the sharp interface between two oth- 
erwise homogeneous porous materials. To this end, we 
begin by resorting to the Continuous Time Random Walk 
(CTRW) formalism, representing the stochastic path of 
each particle as a series of random jumps between spatial 
sites, separated by random waiting times. We denote the 
jump-length probability density function (pdf) by p(s) 
and the waiting time pdf by tp(t), and we assume for the 
sake of simplicity that the two pdfs are uncoupled. The 
pdf ip(t) plays a central role in CTRW, since it defines the 
variability of the velocity spectrum and thus condenses 
the degree of heterogeneity of the traversed medium Q . 
As a first approximation, the two columns can be sepa- 
rately considered as being homogeneous and can thus be 
described by a narrow ip(t) distribution (e.g., an exponen- 
tial pdf), so that a single time scale t (e.g., the mean of 
the distribution) dominates, and particle sojourns at each 
spatial site have on average the same duration. Since the 
two layers have the same porosity, n, we also assume that 
the time scale r is the same in the two layers. 

The flux expression j is based on the assumption 
that the jump lengths are distributed according to a 
Gaussian pdf, g(s) s 
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value of <?(s), //, characterizes the advective contribution 
to the flow, whereas tr, the square root of the variance, 
characterizes the dispersion, a = ^/2ajl. The presence 
of the macroscopic interface imposes a spatial variability 
on the value of ct, but not on the value of [i because the 
porosity in the two sections is the same, so that we can 
make the identification, p(s\x) = g(s\fi, a(x)). This pdf 
reproduces the Fickian flux j AD , and consistently yields 



identical solute arrival times at the outlet upon reversal 
of the flow direction. 

This classical picture of transport needs now to be 
modified to account for the local effects of the interface 
reported in [l[ . Consider a Lagrangian coordinate frame 
moving with velocity v = /u/r. A particle located in 
the fine (coarse) homogeneous material at a sufficiently 
large distance from the interface must have an identical 
probability of jumping in either direction: in this region, 
p(s\x) is symmetrical and entirely characterized in terms 
of v and (either) cv^ (fine layer) or ac (coarse layer). In 
the transition region X, however, a particle has a finite 
probability of starting and ending the jump in two dis- 
tinct layers, and will thus be subject to an asymmetric 
microscopic random force field. For the case of F^C 
flow direction, the probability of jumping forward will be 
larger than the probability of jumping backwards. Such 
a particle will thus locally experience a forward drift, i.e., 
a positively skewed jump length distribution. The same 
reasoning holds, but with opposite signs, for a particle 
crossing the interface in the C— >F direction; such a par- 
ticle will experience a backward drift, i.e., a negatively 
skewed jump length distribution. Skewed Brownian mo- 
tion schemes have been proposed when the forward or 
backward direction of the jumps can be characterized by 
a biased Poisson process (see, e.g., Il9l|) . and for flow 
parallel to a sharp interface (see, e.g., [20(). 

The microscopic nature of the jump length pdf asym- 
metry can be appreciated by observing the behavior of 
the pore-scale Stokes flow velocity field, u, across the 
sharp interface. It is well known from homogenization 
theory that in homogeneous porous media deviations of 
the Stokes velocity from its mean, u' = u — (u), are 
proportional to the permeability k 2l|. Across the in- 



terface of a composite medium we have therefore that 
so that the u' pdf must be skewed. Figure [1] 
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FIG. 3: Breakthrough curves (BTCs) corresponding to three 
values of velocity v = 0.25, 0.5 and 1. Dispersivities are 
of = 0.1 and ac — 0.8, and 9 — 0.35. Monte Carlo simula- 
tion results are displayed as symbols: squares correspond to 
F^C flow conditions, crosses to C— >F. Comparisons with the 
rectified flux model are shown as solid lines. For the sake of 
clarity, the range of Peclet numbers used in these simulations 
is roughly ten times larger than for the experiments in 



shows a simple 2D Stokes flow velocity field finite ele- 
ments simulation [22I ] at the transition between fine and 
coarse disks arrangements: the histogram of the v! com- 
ponent along the flow direction is positively (negatively) 
skewed in the F— >C (C— >F) flow direction. 

A simple, yet realistic way to incorporate this be- 
havior into the functional form of p{s\x G X) is to 
assume that the shape of p(s\x) gradually changes its 
skewness within the region X, in such a way that 
the skewness vanishes when particles are sufficiently 
far from the interface. An expedient pdf that repro- 
duces these features is the skew-normal distribution [16| , 
ff s ) = -L e -(s-v) 2 /^ 2 [ eis -^ /a e ~y 2 / 2 dy, where 9 deter- 

J V ' (T7T J— OO I 1 a 1 

mines the skewness of the pdf [16j. Positive (negative) 
values of 9 correspond to right (left) skewed pdfs, whereas 
for 9 = one recovers the normal (symmetrical) pdf, 
f{s\6 = 0) = g(s). The mean of f(s) equals (s) = p, + £, 
9 [2 

where £ = a , =\ — . Without loss of generality, we 

can characterize the excursion of 9 = 9{x) in the tran- 
sition region by its maximum value 9, positive in the 
F^C direction, and negative in the C^F direction, so 
that p(s\x) = f(s\p,cr(x),0(x)). For a visual represen- 
tation of the evolution of the jump-length skew-normal 
pdfs across the transition region, see Fig. [2] 

The skewness of p(s\x) enters the expression for the 
mean, and thus directly induces a drift velocity correc- 
tion localized at the interface, v' = £/t, which can be 
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FIG. 4: Resident concentration profiles corresponding to 
times t = 0.5, 0.75 and 1, at v = 0.5. Dispersivities are 
olf — 0.01 and ac = 0.08, and 9 = 0.35. Monte Carlo sim- 
ulation results are displayed as symbols: squares correspond 
to F^C flow conditions, crosses to C— >F. Comparisons with 
the rectified flux model are shown as solid lines. 

expressed as 

v> = ±2X9 J \ Q(1/2 L(. - 1) (1) 
V 7r(l + 6> 2 ) V 2> y 1 

where ± signs correspond to the F— >C and C— >F direc- 
tions, respectively, and A = 1 has dimensions of length 
over square root of time. 

We can interpret this rectification flux as arising from 
subtle ratchet-like mechanisms at work in the interface 
between the two layers, X. In other words, small pertur- 
bations in the potential energy of the particles in X do 
not simply average out, but rather induce a succession of 
asymmetric potential wells where the particles can jump 
with a preferential direction. Note that at high veloc- 
ities the particles have enough energy to overcome the 
potential barriers, i.e., the advective part of the solute 
flux dominates over the dispersive part, and the differ- 
ences between the BTCs in the two directions disappear. 
A variety of rectification fluxes have been experimentally 
observed in other systems where nonequilibrium fluctu- 
ations (endogeneous and/or exogeneous) in anisotropic 
media induce a unidirectional bias in the Brownian mo- 
tion of a particle (see, e.g., [l^jEl and refs. therein). 

In Fig. [3J we display one-dimensional Monte Carlo ran- 
dom walk simulations with skew- normal pdfs p(s) = f(s) 
for the breakthrough curves (BTC) measured at the exit 
of the column. Three sets of F^C (squares) and C^F 
curves (crosses) are shown, for increasing values of the 
velocity v. Our simulations qualitatively reproduce the 
salient features experimentally observed in [l| : in partic- 
ular, (i) the C-^F BTC is delayed with respect to the 
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F^C BTC; (m) as v increases, the BTC curves become 
progressively closer to each other (i.e., the delay van- 
ishes), and the standard AD behavior is recovered. 

Monte Carlo estimates of the velocity correction at the 
interface, v', are in good agreement with the predictions 
in Eq. Our model is robust to changes in the func- 
tional form of the p(s\x), as long as the second moment is 
finite and a spatial variation of the skewness is preserved. 

Concentration profiles along the column's longitudi- 
nal direction provide interesting clues about the nature 
of this transport process. The concentration profiles in 
Fig. [4] display significant mass accumulation and sharp 
gradients at the interface, two characteristic signatures 
of our model. A non-invasive concentration profile lab- 
oratory measurement may be thus designed to validate 
the physical mechanism proposed in this work. 

The expression for the solute flux can be written as 
j'(x,t) = j AD (x,t) + v'(x)c(x,t), which, owing to an 
integration by parts argument and neglecting the term 
d y S(y — \), can also be recast as j'(x,t) = j AD (x,t) + 
fo v> (y)dyc(y,t) dy. The equivalence between these two 
flux expressions shows how the interface correction is nei- 
ther purely advective, nor purely dispersive, as also ap- 
parent from our random walk BTCs. We can now numer- 
ically solve the partial differential equation (pde) result- 
ing by inserting the expression for j'(x, t) into the conti- 
nuity equation, and compare the results with the Monte 
Carlo simulations. This is done in Fig. [3] for the BTCs 
and in Fig. 3] for the concentration profiles. The overall 
agreement between the Monte Carlo and pde approaches 
is excellent. From Eq. (fT]), we note that the rectified flux 
depends on the direction of the flow, being positive in 
the F^C direction and negative in the opposite C^F 
direction. The magnitude of the correction flux is small 
when dispersivity contrast is small, and/or the velocity 
v is high. Also, the correction is smaller where the con- 
centration gradient is smaller, i.e., at saturation. Since 
the thickness of the physical interface between the two 
layers strongly depends on the fine details of the Stokes 
velocity field, we expect 9 to increase with the inverse 
of the velocity v, proportionally to the dispersivity and 
permeability contrast. 

When the velocity v decreases below some given 
threshold, diffusion dominates over dispersion: the asym- 
metric transport mechanism can readily be adapted to 
describe this situation. Depending on the specific val- 
ues of the diffusion coefficients for the two homogeneous 
sections, our model predicts that the separation between 
BTCs could even undergo an inversion, i.e., the F^C 
curve could be delayed with respect to the C^F. This is 
a pertinent issue, e.g., for nuclear waste storage in multi- 
layered geological formations, where diffusion is (almost 
everywhere)the dominant transport mechanism. Differ- 
ently from |4j , our model preserves the continuity of the 
resident concentration at the interface. 

Finally, we can extend our model to the case in which 



the two sections are characterized by a tjj(t) pdf with 
slower than exponential decay (small degree of disor- 
der) [l5j ]. Owing to the uncoupling between p(s) and 
i/j(t), the asymmetric transport mechanism at the inter- 
face remains unchanged, and the resulting pde will be 
defined through the time-convolution of j'(x,t) with a 
memory function kernel depending only on ip(t) 
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